Spatiotemporal complexity of a ratio-dependent predator-prey system 
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In this paper, we investigate the emergence of a ratio-dependent predator-prey system with Michaelis-Menten- 
type functional response and reaction-diffusion. We derive the conditions for Hopf, Turing and Wave bifurcation 
on a spatial domain. Furthermore, we present a theoretical analysis of evolutionary processes that involves 
organisms distribution and their interaction of spatially distributed population with local diffusion. The results 
of numerical simulations reveal that the typical dynamics of population density variation is the formation of 
isolated groups, i.e., stripelike or spotted or coexistence of both. Our study shows that the spatially extended 
model has not only more complex dynamic patterns in the space, but also chaos and spiral waves. It may help 
us better understand the dynamics of an aquatic community in a real marine environment. 

PACS numbers: 87.23.Cc, 89.75.Kd, 89.75.Fb, 47.54. -r 



INTRODUCTION 

Ecological systems are characterized by the interaction be- 
tween species and their natural environment m. Such inter- 
action may occur over a wide range of spatial and temporal 
scales J2, The study of complex population dynamics 
is nearly as old as population ecology. In the 1920s, Lotka 
and Volterra independently developed a simple model of in- 
teracting species that still bears their joint names. This was 
a nearly linear model, but the predator-prey version displayed 
neutrally stable cycles 0,|5|. From then on, the dynamic re- 
lationship between predators and their prey has long been and 
will continue to be one of dominant themes in both ecology 
and mathematical ecology due to its universal existence and 
importance JEHU- 

Predator-prey models follow two general principles: one is 
that population dynamics can be decomposed into birth and 
death processes; the other is the conservation of mass princi- 
ple, stating that predators can grow only as a function of what 
they have eaten J<J. With these two principles we can write 
the canonical form of a predator-prey system as 



x(t) = xg{x) - f{x,y)y - fi x (x)x, 
y{t) = if{x,y)y - v v {y)y- 



(l) 



where g(x) is the per capita prey growth rate in the absence 
of the predator, \i x and fj, y are natural mortalities of prey and 
predator respectively, /(x, y) is the functional response. And 
7/(x, y) is the per capita production of predator due to preda- 
tion, which is often called the numerical response. The func- 
tional response plays a main role in system (Q3: the knowledge 
of this function determines the dynamics of the whole system 
and the transfer of the biomass in the predation because it is 
proportional to the numerical response. Usually one consid- 
ers consumption to be the major death cause for the prey. In 
this case \i x (x) can be neglected and set to (as long as the 
predator exists) Jit] . 

In population dynamics, a functional response of the preda- 
tor to the prey density refers to the change in the density of 



prey attached per unit time per predator as the prey density 
changes ifioll . In general, functional response can be classified 
as: (i) prey dependent, when prey density alone determines the 
response, i.e., f(x, y) = p(x); (ii) predator dependent, when 
both predator and prey populations affect the response. Partic- 
ularly, when /(x, y) = p(^), we call model (Q~|) strictly ratio- 
dependent; and (iii) multi-species dependent, when species 
other than the focal predator and its prey species influence the 
functional response ll ill . Differing from the prey-dependent 
predator-prey models, the ratio-dependent predator-prey sys- 
tems have two principal predictions: (a) equilibrium abun- 
dances are positively correlated along a gradient of enrich- 
ment and (b) the "paradox of enrichment" either completely 
disappears or enrichment is linked to stability in a more com- 
plex way II121 1 1 3f1 - The ratio-dependent predator-prey model 



has been studied by several researchers recently and very rich 
dynamics have been observed if/l, lil [12 , 14, 151. 



On the other hand, pattern formation in nonlinear complex 
systems is one of the central problems of the natural, social, 
and technological sciences 111 611 . In particular, starting with 
the pioneering work of Segel and Jackson 01711 . spatial pat- 
terns and aggregated population distributions are common in 
nature and in a variety of spatio-temporal models with local 
ecological interactions ifU flal . Promulgated by the theoretical 
paper of Turing lfl9Tl . the field of research on pattern forma- 
tion modeled by reaction-diffusion systems, which provides 
a general theoretical framework for describing pattern forma- 
tion in systems from many diverse disciplines including (but 

" 23, EI 



not limited to) biolo 
istry ill |27l EH 



25], chem- 
3511 . and so 



lia 1201 [2JJ, [22L123L |2; 
30j|3l|], physics Il32l.l33l. l3 

on, seems to be a new increasingly interesting area, particu- 
larly during the last decade. 

In Ref. J2|, Neuhauser surveys some current work on spatial 
mathematical models in ecology. Much of this work consists 
of building spatial dimensions into existing classical models, 
such as the Lotka- Volterra model that describes competition 
between species. But the research on the spatial pattern of 
ratio-dependent predator-prey models seems to be rare. 
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STABILITY AND BIFURCATION ANALYSIS 

In this paper, we mainly focus on the ratio-dependent 
predator-prey system with Michaelis-Menten-type(or 
Michaelis-Menten-Holling) functional response: 



M. _ /-i _ JVn m _ aN/P 
dt — 1 \ L Kl 1 ^ 1+ahN/P 



= r(l-%)N 



aN 



9P =y _ 
dt I P+ahN 



P+ahN 

aN P-fiP + D 2 V 2 P. 

2 



P + D{V 2 N 
P + D{V 2 N, 



(2) 



I V(jV,P) g [0,^] 2 \(0,0). 



where N, P stand for prey and predator density, respectively. 
Di, D 2 are their respective diffusion coefficients, V 2 = + 
g|a . All parameters are positive constants, r stands for maxi- 
mal growth rate of the prey, 7 conversion efficiency, /j, preda- 
tor death rate, K carrying capacity, a capture rate and h han- 
dling time. 

Note that 1+ a < ^J^ P is strictly correct only for P > 0. In 
the case of P = and N > we can define f(N, 0) := h 



(the limit of fix) for x 
Let 



00). 



N 



Q 



p 



ahN 
~ ( K > 

h[i cf 

— 1 — 



ahN 



R = ^, 

7 ' 



t = 2* 



(3) 



For simplicity we will not write the hat (") in the rest of this 
paper. And in these new variables, from ©-(O, we arrive at 
the following equations containing dimensionless quantities: 



dN _ p>(i _ N_\at _ _ 
dt ^H 1 s> P+SN 



SN 



P + D^N, 



dP _ SN 

at p+sn 



P-QP + D 2 V 2 P. 



(4) 



More details about the choice of dimensionless variables in 
the system (f2]i as well as possible implications can be found 
in H. 

The dimensionless model (Eq.|4]i has only three parameters: 
R, which controls the growth rate of prey; Q, which controls 
the death rate of the predator; and S, which measures captur- 
ing rate. 

The first step in analyzing the model is to determine the 
behavior of the non-spatial model obtained by setting space 
derivatives equal to zero. The non-spatial model has at most 
three equilibria (stationary states), which correspond to spa- 
tially homogeneous equilibria of the full model (Eq. |4), in 
the positive quadrant: (0, 0) (total extinct), (5, 0) (extinct of 
the predator) and a nontrivial stationary state (n* ,p*) (coex- 
istence to prey and predator), where 



» _ s(BMQ-i)S) 

11 ~ R ' 

» _ S(l-Q) « _ S 2 (R-S+QS){1-Q) 
p Q 11 RQ 



(5) 



Easy to know that n* is positive for all S < j^q, which 



implies Q < 1 and therefore ensures the positivity of p* IU4I1 . 



To perform a linear stability analysis, we linearize the 
dynamic system (O around the spatially homogenous fixed 
point (0 for small space- and time-dependent fluctuations and 
expand them in Fourier space 



N(x,t) - n*e xt e lk - s , 
P(x,t) ~p*e xt e zl3 . 
and obtain the characteristic equation 

\A - k 2 D - \I\ =0, 

where 

°=( Di °v 

V D 2 J 

and A is given by 

/ d N f d P f \ ( }n fp 



(6) 



(7) 



(8) 



9n9 d P g 



(9) 



(n*,p*) 



9n gp 



where the elements are the partial derivatives of the reaction 
kinetics evaluated at the stationary state (n*,p*). Now Eq. (0 
can be solved, yielding the so called characteristic polynomial 
of the original problem (Eq. |4]i 



X 2 -tr k X + A k = 0, 



(10) 



where 



tr k =f N +gp-k 2 (D 1 +D 2 ) =tro-fc 2 0Di+A2), (11) 

A fc = f N gp - fpg N - k 2 {f N D 2 + gpDt) + k^D x D 2 
= A - k 2 {f N D 2 + g P D x ) + k i D 1 D 2l 



The roots of Eq. 10 yield the dispersion relation 
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(12) 



(13) 



It's well known that reaction-diffusion systems have led 
to the characterization of three basic types of symmetry- 
breaking bifurcations responsible for the emergence of spa- 
tio temporal patterns. The space-independent Hopf bifurcation 
breaks the temporal symmetry of a system and gives rise to os- 
cillations that are uniform in space and periodic in time. The 
(stationary) Turing bifurcation breaks spatial symmetry, lead- 
ing to the formation of patterns that are stationary in time and 
oscillatory in space. The wave (oscillatory Turing or finite- 
wavelength Hopf) bifurcation breaks both spatial and tempo- 
ral symmetry, generating patterns that are oscillatory in space 
and time Q. 

The Hopf bifurcation occurs when 

Jm(A(fc)) + 0, Re(X(k)) = at k = 0. (14) 
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then we can get the critical value of Hopf bifurcation parame- 
ter S equals 



5fl 



R + Q-Q 2 
1-Q 2 



(15) 



At the Hopf bifurcation threshold, the temporal symmetry of 
the system is broken and gives rise to uniform oscillations in 
space and periodic oscillations in time with the frequency 



cj h = Im(\(k)) = ^Ao = y/Q(Q-l)(R-S + QS), 
the corresponding wavelength is 



A H = 



2tt 



2 k 



(16) 



y/Q(Q-l)(R-S + QS) 
The Turing bifurcation occurs when 
Im(A(fc)) = 0, Re(X(k)) = at k = k T ± 0. (17) 
the critical value of bifurcation parameter S equals 

D x D 2 k T + {D 2 R + D X Q(1 - Q))k 2 T + RQ{1 - Q) 



where 



Q3 _ ( k 2 T D 2 + 2)Q 2 + Q + k 2 T D 2 



(18) 



D X D 2 



and at the Turing threshold, the spatial symmetry of the sys- 
tem is broken and the patterns are stationary in time and os- 
cillatory in space with the wavelength 



at — - — ■ 



(19) 



And the Wave bifurcation occurs when 

Jm(A(fc)) ± 0, Re(X(k)) = at k = k w ^ 0. (20) 
the critical value of Wave bifurcation parameter S equals 

k 2 w {D l +D 2 )+R + Q-Q 2 



Sw 



1 



(21) 



where 



- ) {{D l -D 2 ) 2 Q i -2{D 2 +D 2 )Q 3 



k 2 = £ 

w 2DKQ+ 



-(D\ + Dl+ioD^ - 4D%R)Q 2 - AD 1 D 2 Q + ±D\R 



1/2 



Easy to know that, at the Wave threshold, both spatial and tem- 
poral symmetries are broken and the patterns are oscillatory in 
space and time with the wavelength 



A - 2n 

AW — -, — • 



(22) 



Linear stability analysis yields the bifurcation diagram with 
R = 0.5, Q = 0.6, D 2 = 0.2 shown in Fig. Q] 



TlUllLg 
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FIG. 1: (Color online) Bifurcation diagram for the system (2) with 
R = 0.5, Q = 0.6, D 2 = 0.2. Hopf bifurcation line: Sh = §f; 
Turing bifurcation line: St — — 1 

bifurcation line: Sw = (ff-Di TE) K w 
Turing-Hopf bifurcation point: (0.08864, 1.15625) and Turing- Wave 
bifurcation point: (0.11655, 1.21110) 



56 +H Di with = 5.48571; Wave 
+ § with k w = 0.334. 



The Hopf bifurcation line, the Wave bifurcation line and 
the Turing bifurcation line intersect at two codimension-2 bi- 
furcation points, the Turing-Hopf bifurcation point and the 
Turing-Wave bifurcation point. The bifurcation lines sepa- 
rate the parametric space into six distinct domains. In domain 
I, located below all three bifurcation lines, the steady state is 
the only stable solution of the system. Domain II is region of 
pure Turing instabilities, and domain III is pure Hopf instabil- 
ities. In domain IV, both Hopf and Turing instabilities occur, 
and in domain V, the Wave and Hopf modes arise. When the 
parameters correspond to domain VI, which is located above 
all three bifurcation lines, all three instabilities occur. Fig. [2] 
shows the dispersion relations of unstable Hopf mode, transi- 
tion of Turing and Wave modes from stable to unstable. It's 
easy to know that all three bifurcations are supercritical. 



SPATIOTEMPORAL PATTERN ANALYSIS 

We have performed extensive numerical simulations of the 
spatially extended model (01 in two-dimensional space, and 
the qualitative results are shown here. All our numerical sim- 
ulations employ the periodic Neumann (zero-flux) boundary 
conditions with a system size of 200 x 200 space units and 
R = 0.5, Q = 0.6, Di = 0.02, D 2 = 0.2. The spatiotem- 
poral dynamics of a diffusion-reaction system depends on the 
choice of initial conditions, which some authors have consid- 
ered in connection with t he p roblem of biological invasion in 
a few papers |T(| 36[ 37, 38], where the initial conditions are 
naturally described by finite functions and the dynamics of the 
community mainly consists of a variety of diffusive popula- 
tional fronts. In general, there are two initial conditions used 
for analysis of the spatial extended systems. One is random 
spatial distribution of the species, which seems to be more 
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FIG. 2: (Color online) Dispersion relations showing unstable Hopf 
mode, transition of Turing and Wave modes from stable to unsta- 
ble, e.g., as Di decreased. Parameters: S — 1.2, R = 0.5, Q = 
0.6, D 2 = 0.2 and (1) £>i=0.15; (2) Di=0.12; (3) L>i=0.10; (4) 
£>i=0.07; (5) Di=0.04; (6) £>i=0.02. 



general from the biological point of view (cf. Fig. [3jA) and 
EAX |2lA)); The other is a special choice, i.e., taking the 
species community in a horizontal layer as decreasing gradu- 
ally and the vertical distribution of species homogeneous(cf. 
Fig. IDA)). In this section we choose the former, and the lat- 
ter in the next section (IV. Discussion). The equations (|4]i 
are solved numerically in two-dimensional space using a fi- 
nite difference approximation for the spatial derivatives and 
an explicit Euler method for the time integration with a time 
stepsize of At = 0.01 and space stepsize Ah = 0.25. 

From the analysis of section II and phase-transition bifurca- 
tion diagram (cf. Fig.[T|), the results of computer simulations 
show that the type of the system dynamics is determined by 
the values of S and D\, We run the simulations until they 
reach a stationary state or until they show a behavior that does 
not seem to change its characteristics anymore. For different 
sets of parameters, the features of the spatial patterns become 
essentially different if S exceeds a critical value St, St and 
Sw respectively, they depend on D\. 

Fig.[3]shows the evolution of the spatial pattern of prey at 0, 
5000 and 45000 iterations, with random small perturbation of 
the stationary solution n* andp* of the spatially homogeneous 
systems when £ is less than the Turing bifurcation threshold 
St- In this case, one can see that for the system the ran- 
dom initial distribution lead to the formation of a strongly ir- 
regular transient pattern in the domain. After the irregular 
pattern form (cf. Fig.[2B)) it grows slightly and "jumps" al- 
ternately for a certain time, and finally the chaos spiral pat- 
terns prevail over the whole domain, and the dynamics of the 
system does not undergo any further changes (cf. Fig. OQ 
and the addition movie for Fig. [3]). 

Figures H] and [5] show spontaneous formation of short 
stripelike and spotted spatial patterns emerge and coexist sta- 
bly when the bifurcation parameter S < St (Fig. HJi and 
St < S < Sh (Fig. |5]l. From the snapshots or movies, one 




FIG. 3: (Color online) Snapshots of contour pictures of the time evo- 
lution of the prey at different instants with S = 0.6 < St- (A) 
iteration; (B) 5000 iterations; (C) 45000 iterations. [Additional 
movie format available from the author] 



can see that the stripelike spatial patterns arise from the ran- 
dom initial conditions. After the stripelike patterns form (cf. 
Fig. HtB)) they grow steadily with time until they reach cer- 
tain width — armlength, and the spatial patterns become dis- 
tinct. Finally, the stripelike spatial patterns prevail the whole 
domain (cf. Fig. HJC) and FigfS). Comparing the Fig. EfC) 
with Fig.|5j we find that the parameter 5 is closer to St, and 
the stripelike spatial patterns are more distinct. Here we omit 
the pre-image of Fig. [5] as they are similar to the Fig.|4|A) and 
(B). It is easy to see that the stationary patterns are essentially 
different from the previous case(cf. Fig[3]). 




FIG. 4: (Color online) Snapshots of contour pictures of the time evo- 
lution of the prey at different instants with S = 0.9 < St- (A) 
iteration; (B) 5000 iterations; (C) 45000 iterations. [Additional 
movie format available from the author] 




FIG. 5: (Color online) Snapshots of contour pictures of the time evo- 
lution of the prey at different instants with St < S = 1.1 < Sh 
(45000 iterations). [Additional movie format available from the au- 
thor] 

When S H < S = 1.16 < S w , we find that the spotted 
patterns and the stripelike patterns coexist in the spatially ex- 
tended model (cf. Fig.|6]l. 

Fig.|7]shows snapshots of prey spatial pattern at 0, 5000 and 
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the original images. In Fig.|8jmid-hand column), short wave- 
length, represented by a large circle, corresponds to Turing 
structures; longer wavelength, represented by a small white 
circle, corresponds to traveling and/or standing waves. This 
technique can be particularly appropriate for characterizing 
quasi-ordered arrays for Fig.|7jC). 



DISCUSSION AND CONCLUSION 



FIG. 6: (Color online) Coexistence of stationary spotted patterns and 
stripelike patterns of the prey for long time run with Sh < S = 
1.16 < Sw- [Additional movie format available from the author] 



45000 iterations for the parameter S = 1.2 > Sw- Although 
the dynamics of the system starts from the same initial condi- 
tion as previous cases, there is an essential difference for the 
spatially extended model (Eq. @}. Form Fig. [7] one can see 
that the regular spotted patterns prevails over the whole do- 
main at last, and the dynamics of the system does not undergo 
any further changes (cf. the additional movie for Pig - EJ . 




FIG. 7: (Color online) Snapshots of contour pictures of the time evo- 
lution of the prey at different instants with S = 1.2 > Sw- (A) 
iteration; (B) 5000 iterations; (C) 45000 iterations. [Additional 
movie format available from the author] 

On the other hand, discrete Fourier transform is a basic 
mathematical tool used to decompose a signal or image into 
different periodic c omp onents. It has been widely used for 
11E 



the spatial patterns [3 



41H . We have also performed nu- 



merical investigations into two-dimensional space by Fourier 
spectra. The numerical computation of the Fourier transform 
is done by the well-established two-dimensional Fast Fourier 
Transform (FFT2) algorithm 14211 . Spatial Fourier transform 
of the stripelike and spotted patterns in Figures SJc), [5] and 
Ulc) are shown as Fig. [8] And digital or digitized transmission 
electron micrographs (TEMs) are analyzed by using Matlab 
(Ver.7.0). 

From Fig.(8j we find that Fig.[4|C) an d Fig.[5]have the same 
spatial frequency in the length of the space unit and presence 
of one mode with different wave numbers. On the contrary, 
Fig- EtC) has two modes with different wave numbers. The 
spatial frequency and direction of any component in the power 
spectrum are given in the length and direction, respectively, 
of a vector from the origin to the point on the circle. The 
magnitude is depicted by a gray scale or color scale, but the 
units are dimensionless values related to the total darkness of 



The numerical results correspond perfectly to our theoret- 
ical findings that there are a range of parameters in S — D\ 
plane where the different spatial patterns emerge (cf. Fig.Q]). 
Fig.Q]and the results of simulations present that the chaos pat- 
terns will persist in the spatially extended model (Eq.[4]i when 
the parameters are in the domain I. The boundary of this do- 
main can be computed numerically and is shown as the blue 
line "Turing" in Fig.Q] The stationary state of stripelike pat- 
terns exists when the parameters are in the domain II, where 
its boundary can also be computed numerically and is shown 
as the black line "Hopf" in Fig. Q] The periodic spotted pat- 
terns appear in the domain VI, where the boundary can be 
computed numerically by the wave bifurcation and is shown 
as the red line "Wave" in Fig.Q] Moreover, there is transverse 
domain IV (cf. Fig. [TJ in the system between the stripelike 
patterns and spotted patterns, where the spotted patterns and 
the stripelike patterns coexist(cf. Fig.|6]l. 

Do the stationary patterns arise dependent on the initial con- 
ditions? We test the different initial conditions for the spatially 
extended system, but the final spatial patterns are the same 
in qualitative. In those figures we find that the spatial chaos 
patterns come from the destruction of the spirals, when we 
choose the special initial condition (cf. Fig.|9]l in the domain 
I. This phenomenon coheres with the results of the study in 
Refs. dlH. 

We have presented a theoretical analysis of evolutionary 
processes that involves organisms distribution and their inter- 
action of spatially distributed population with local diffusion. 
Our analysis and numerical simulations reveal that the typical 
dynamics of population density variation is the formation of 
isolated groups (stripelike or spotted or coexistence of both). 
This process depends on several parameters, including S, D\ 
and Z?2- The field meaning of our results may be found in the 
dynamics of an aquatic community which is affected by the 
existence of relatively stable mesoscale inhomogeneity in the 
field of ecologically significant factors such as water temper- 
ature, salinity and biogen concentration. 

In Ref. 11611 . the authors explained the field meaning by us- 
ing aquatic community in the ocean (cf. Fig. [5). Our study 
shows that the spatially extended model (Eq. [4]i has not only 
more complex dynamic patterns in the space, but also chaos 
patterns and spiral waves, so it may help us better understand 
the dynamics of an aquatic community in a real marine envi- 
ronment. It is also important to distinguish between "intrin- 
sic" patterns, i.e., patterns arising due to trophic interactions 
like those considered above, and "forced" patterns induced by 
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FIG. 8: (Color online) Patterns (left-hand column), it spatial Fourier transformation (mid-hand column), and radial average of the power 
spectrum (right-hand column). Up-row S = 0.9; Mid-row S = 1.1; Low-row S = 1.2. 




FIG. 9: (Color online) Snapshots of contour pictures of the time evo- 
lution of the prey at different instants with the special initial condition 
and 5" = 0.6 < St. (A) iteration; (B) 1000 iterations; (C) 400000 
iterations. [Additional movie format available from the author] 



ratio-dependent predator-prey model (Eq. |4]i also represents 
rich spatial dynamics, such as chaos spiral patterns, stripelike 
patterns, spotted patterns, coexistence of both stripelike and 
spotted patterns, etc. It will be useful for studying the dy- 
namic complexity of ecosystems. 

This work was supported by the National Natural Science 
Foundation of China under Grant No. 10471040 and the 
Natural Science Foundation of Shan'xi Province Grant No. 
2006011009. 



the inhomogeneity of the environment. The physical nature 
of the environmental heterogeneity, and thus the value of the 
dispersion of varying quantities and typical times and lengths, 
can be essentially different in different cases. Neuhauser and 
Pacala |44J| formulated the Lotka-Volterra model as a spatial 
model. They found the striking result that the coexistence of 
patterns is actually harder to get in the spatial model than in 
the non-spatial one. One reason can be traced to how local in- 
teractions between individual members of the species are rep- 
resented in the model. In this thesis, our results show that the 
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